R Notebook

The goal here is to use conventional alpha diversity metrics to see how Chao1 richness, shannon diversity and evenness change across samples and to compare those to the values seen using breakaway in the AlphaDiversity.Rmd file

Setup

Run AlphaDiversity in scratchnotebooks That file calculates richness in breakawy which I will combine here

#source(here::here("RScripts", "InitialProcessing_3.R"))
source(here::here("RLibraries", "ChesapeakePersonalLibrary.R"))

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.4.0      ✔ purrr   0.3.4 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.3      ✔ forcats 0.5.2 ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

ksource(here::here("ActiveNotebooks", "BreakawayAlphaDiversity.Rmd"))

processing file: /home/jacob/Projects/ChesapeakeMainstemAnalysis_ToShare/ActiveNotebooks/BreakawayAlphaDiversity.Rmd

  |                                                                                                                                                                            
  |                                                                                                                                                                      |   0%
  |                                                                                                                                                                            
  |.......                                                                                                                                                               |   4%
  |                                                                                                                                                                            
  |..............                                                                                                                                                        |   9%
  |                                                                                                                                                                            
  |......................                                                                                                                                                |  13%
  |                                                                                                                                                                            
  |.............................                                                                                                                                         |  17%
  |                                                                                                                                                                            
  |....................................                                                                                                                                  |  22%
  |                                                                                                                                                                            
  |...........................................                                                                                                                           |  26%
  |                                                                                                                                                                            
  |...................................................                                                                                                                   |  30%
  |                                                                                                                                                                            
  |..........................................................                                                                                                            |  35%
  |                                                                                                                                                                            
  |.................................................................                                                                                                     |  39%
  |                                                                                                                                                                            
  |........................................................................                                                                                              |  43%
  |                                                                                                                                                                            
  |...............................................................................                                                                                       |  48%
  |                                                                                                                                                                            
  |.......................................................................................                                                                               |  52%
  |                                                                                                                                                                            
  |..............................................................................................                                                                        |  57%
  |                                                                                                                                                                            
  |.....................................................................................................                                                                 |  61%
  |                                                                                                                                                                            
  |............................................................................................................                                                          |  65%
  |                                                                                                                                                                            
  |...................................................................................................................                                                   |  70%
  |                                                                                                                                                                            
  |...........................................................................................................................                                           |  74%
  |                                                                                                                                                                            
  |..................................................................................................................................                                    |  78%
  |                                                                                                                                                                            
  |.........................................................................................................................................                             |  83%
  |                                                                                                                                                                            
  |................................................................................................................................................                      |  87%
  |                                                                                                                                                                            
  |........................................................................................................................................................              |  91%
  |                                                                                                                                                                            
  |...............................................................................................................................................................       |  96%
  |                                                                                                                                                                            
  |......................................................................................................................................................................| 100%

output file: /tmp/RtmphhTqE3/file42a957179148

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio

Attaching package: ‘flextable’

The following object is masked from ‘package:purrr’:

    compose


Attaching package: ‘ftExtra’

The following object is masked from ‘package:flextable’:

    separate_header

Warning: Assuming taxa are rows

library(vegan)

Loading required package: permute
Loading required package: lattice
This is vegan 2.6-3

library(cowplot)
library(flextable)
library(ftExtra)

This file is dedicated to conventional, non div-net/breakaway stats, since breakaway seems to choke on this data.

Reshape back into an ASV matrix, but this time correcting for total abundance

raDf <- nonSpikes_Remake %>% pivot_wider(id_cols = ID, names_from = ASV, values_from = RA, values_fill = 0)
raMat <- raDf %>% column_to_rownames("ID")

raMat1 <- raMat %>% as.matrix()

countMat <-  nonSpikes_Remake %>%
  pivot_wider(id_cols = ID, names_from = ASV, values_from = reads, values_fill = 0) %>%
  column_to_rownames("ID") %>% as.matrix()

seqDep <- countMat %>% apply(1, sum)
min(seqDep)

[1] 852

sampleRichness <- rarefy(countMat, min(seqDep))

rarefy everything to the minimum depth (852)

countRare <- rrarefy(countMat, min(seqDep))

Gamma diversity

specpool(countRare)

Doesn’t finish

#specpool(countMat)

Calculate diversity indeces

All richness estimates

richnessRare <- estimateR(countRare)

Shannon diversity

shan <- diversity(countRare)
shan

 3-1-B-0-2  3-1-B-1-2  3-1-B-180   3-1-B-20    3-1-B-5  3-1-B-500   3-1-B-53  3-1-S-0-2  3-1-S-1-2  3-1-S-180   3-1-S-20    3-1-S-5  3-2-B-0-2  3-2-B-1-2  3-2-B-180   3-2-B-20 
  4.377867   5.122635   4.697075   5.919422   5.091001   3.868273   5.527489   4.512680   4.845492   4.782770   5.290467   4.835734   4.308318   4.664723   4.653949   5.022267 
   3-2-B-5  3-2-B-500   3-2-B-53  3-2-S-0-2  3-2-S-1-2  3-2-S-180   3-2-S-20    3-2-S-5  3-2-S-500   3-2-S-53  3-3-B-0-2  3-3-B-1-2  3-3-B-180   3-3-B-20    3-3-B-5  3-3-B-500 
  5.364248   4.956451   4.844607   3.744475   5.000444   4.749852   4.918364   5.142239   4.981104   4.204087   4.331178   4.789362   3.432382   5.736709   5.205365   5.265310 
  3-3-B-53  3-3-S-180   3-3-S-20  3-3-S-500   3-3-S-53  4-3-B-0-2  4-3-B-1-2  4-3-B-180   4-3-B-20    4-3-B-5  4-3-B-500   4-3-B-53  4-3-O-1-2  4-3-O-180    4-3-O-5  4-3-O-500 
  5.492566   4.997043   4.849876   4.869032   4.255819   4.360230   4.996300   4.423317   4.437111   4.771900   4.308965   4.327765   4.891794   4.683325   5.170850   4.084572 
  4-3-O-53  4-3-S-0-2  4-3-S-180   4-3-S-20  4-3-S-500   4-3-S-53  5-1-S-1-2  5-1-S-180   5-1-S-20    5-1-S-5  5-1-S-500   5-1-S-53  5-5-B-0-2  5-5-B-180    5-5-B-5  5-5-B-500 
  4.491632   2.805419   4.592134   4.741784   4.721501   4.533240   4.428191   4.456642   4.070497   3.924195   4.118489   4.097255   4.642190   5.118735   5.533426   5.073586 
  5-5-B-53  5-5-S-180    5-5-S-5  5-5-S-500   5-5-S-53 C_5P1B_0P2 C_5P1B_180 C_5P1B_1P2  C_5P1B_20 C_5P1B_500  C_5P1B_53 
  4.951162   4.296940   4.977480   4.909048   4.280855   4.186882   4.900147   4.860817   5.399822   4.662039   5.209531 

Evenness

pielouJ <- shan/richnessRare["S.chao1",]
pielouJ

  3-1-B-0-2   3-1-B-1-2   3-1-B-180    3-1-B-20     3-1-B-5   3-1-B-500    3-1-B-53   3-1-S-0-2   3-1-S-1-2   3-1-S-180    3-1-S-20     3-1-S-5   3-2-B-0-2   3-2-B-1-2 
0.011467055 0.007103701 0.012255238 0.003349956 0.008110994 0.041006425 0.007615813 0.011728862 0.009093947 0.008972655 0.008846283 0.007282732 0.012661180 0.009597579 
  3-2-B-180    3-2-B-20     3-2-B-5   3-2-B-500    3-2-B-53   3-2-S-0-2   3-2-S-1-2   3-2-S-180    3-2-S-20     3-2-S-5   3-2-S-500    3-2-S-53   3-3-B-0-2   3-3-B-1-2 
0.011329508 0.005201462 0.008190891 0.009811517 0.010195269 0.021596351 0.008966390 0.010189900 0.008947565 0.008349931 0.009991783 0.014580123 0.012148164 0.008882062 
  3-3-B-180    3-3-B-20     3-3-B-5   3-3-B-500    3-3-B-53   3-3-S-180    3-3-S-20   3-3-S-500    3-3-S-53   4-3-B-0-2   4-3-B-1-2   4-3-B-180    4-3-B-20     4-3-B-5 
0.063562627 0.005282231 0.009346957 0.006143886 0.005836946 0.010690993 0.008613137 0.009301306 0.010889362 0.012524899 0.007823094 0.011012647 0.008667124 0.009030647 
  4-3-B-500    4-3-B-53   4-3-O-1-2   4-3-O-180     4-3-O-5   4-3-O-500    4-3-O-53   4-3-S-0-2   4-3-S-180    4-3-S-20   4-3-S-500    4-3-S-53   5-1-S-1-2   5-1-S-180 
0.005918521 0.010062620 0.009452742 0.014296051 0.006724997 0.014380807 0.011650840 0.124685269 0.011358859 0.008391334 0.010995971 0.011373619 0.011834659 0.012159030 
   5-1-S-20     5-1-S-5   5-1-S-500    5-1-S-53   5-5-B-0-2   5-5-B-180     5-5-B-5   5-5-B-500    5-5-B-53   5-5-S-180     5-5-S-5   5-5-S-500    5-5-S-53  C_5P1B_0P2 
0.014452214 0.011299615 0.014013080 0.013216952 0.009980342 0.008042424 0.005883545 0.010729090 0.008357545 0.011428033 0.011626609 0.008644633 0.015737093 0.008575283 
 C_5P1B_180  C_5P1B_1P2   C_5P1B_20  C_5P1B_500   C_5P1B_53 
0.007366125 0.006551531 0.004008265 0.007772224 0.005160141 

Combine diversity data

diversityData <- sampleData %>% 
  left_join(richnessRare %>% t() %>% as.data.frame() %>% rownames_to_column("ID"), by = "ID") %>%
  left_join(shan %>% enframe(name = "ID", value = "shannonH"), by = "ID") %>%
  left_join(pielouJ %>% enframe(name = "ID", value = "pielouJ"), by = "ID") %>%
  arrange(Size_Class)

Generate plots of diversity estimates

Parameters for all plots

plotSpecs <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

Observed species counts, on rarefied data

plotObs <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.obs, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +ylab("Observed ASVs (Rarefied)")#+ scale_y_log10()
plotObs

Estemated Chao1 Richness

plotChao1 <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = S.chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = S.chao1 -2 * se.chao1, ymax = S.chao1 + 2* se.chao1), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Chao1)")
plotChao1

Shannon diversity

plotShan <- diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = shannonH, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  ylab("Diversity (Shannon H)") +
  lims(y = c(2.5, 6))
plotShan

Evenness

plotPielou <- diversityData %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +scale_y_log10() +ylab("Evenness (PielouJ)")
plotPielou

All plots together

plotAlpha <- plot_grid(plotObs, plotChao1, plotShan, plotPielou, nrow = 1, labels = LETTERS)
plotAlpha

ggsave(here::here("Figures", "ConventionalAlpha.png"), plotAlpha, width = 11, height = 4)

Do we see trends with lat and size?

Observed Species

Rarefied observed species numbers

obsMod <- lm(S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(obsMod)

Call:
lm(formula = S.obs ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
     Min       1Q   Median       3Q      Max 
-219.518  -35.649   -1.039   49.226  201.182 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)          248133.486 132566.735   1.872   0.0655 .  
log(Size_Class)          33.680      8.032   4.193 8.02e-05 ***
I(log(Size_Class)^2)     -6.473      1.494  -4.332 4.92e-05 ***
lat                  -12915.027   6907.881  -1.870   0.0658 .  
I(lat^2)                168.133     89.927   1.870   0.0658 .  
depth                     3.905      3.181   1.227   0.2238    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 75.51 on 69 degrees of freedom
Multiple R-squared:  0.2656,    Adjusted R-squared:  0.2124 
F-statistic: 4.991 on 5 and 69 DF,  p-value: 0.0005849

Richness

Rarified chao1 estimates

chao1Mod <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(chao1Mod)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-521.02 -131.00  -37.98  109.00 1052.27 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)          795248.56  436480.86   1.822 0.072796 .  
log(Size_Class)          89.88      26.45   3.399 0.001128 ** 
I(log(Size_Class)^2)    -17.77       4.92  -3.611 0.000574 ***
lat                  -41431.39   22744.45  -1.822 0.072850 .  
I(lat^2)                539.62     296.09   1.822 0.072715 .  
depth                    16.95      10.47   1.618 0.110139    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 248.6 on 69 degrees of freedom
Multiple R-squared:  0.1986,    Adjusted R-squared:  0.1405 
F-statistic: 3.419 on 5 and 69 DF,  p-value: 0.008097

As above but without latitude and depth

chao1ModSimple <- lm(S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), data = diversityData)
summary(chao1ModSimple)

Call:
lm(formula = S.chao1 ~ log(Size_Class) + I(log(Size_Class)^2), 
    data = diversityData)

Residuals:
    Min      1Q  Median      3Q     Max 
-468.58 -139.15  -22.35  112.75 1118.50 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           539.262     44.634  12.082  < 2e-16 ***
log(Size_Class)        90.573     26.507   3.417 0.001045 ** 
I(log(Size_Class)^2)  -18.060      4.931  -3.663 0.000473 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 249.5 on 72 degrees of freedom
Multiple R-squared:  0.1577,    Adjusted R-squared:  0.1343 
F-statistic:  6.74 on 2 and 72 DF,  p-value: 0.002074

Shannon Diversity

shanMod <- lm(shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)

summary(shanMod)

Call:
lm(formula = shannonH ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.32735 -0.16629  0.02731  0.30484  0.76295 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.419e+03  7.933e+02   1.788   0.0782 .  
log(Size_Class)       2.055e-01  4.807e-02   4.276 6.00e-05 ***
I(log(Size_Class)^2) -3.699e-02  8.943e-03  -4.136 9.79e-05 ***
lat                  -7.362e+01  4.134e+01  -1.781   0.0793 .  
I(lat^2)              9.580e-01  5.382e-01   1.780   0.0795 .  
depth                 1.321e-02  1.904e-02   0.694   0.4899    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4519 on 69 degrees of freedom
Multiple R-squared:  0.2905,    Adjusted R-squared:  0.239 
F-statistic: 5.649 on 5 and 69 DF,  p-value: 0.0002011

Evenness

pielouMod <- lm(pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityData)
summary(pielouMod)

Call:
lm(formula = pielouJ ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityData)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.014081 -0.004516 -0.001758  0.000707  0.099175 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -2.476e+01  2.619e+01  -0.946   0.3476  
log(Size_Class)      -4.146e-03  1.587e-03  -2.613   0.0110 *
I(log(Size_Class)^2)  7.087e-04  2.952e-04   2.401   0.0191 *
lat                   1.289e+00  1.365e+00   0.945   0.3480  
I(lat^2)             -1.676e-02  1.776e-02  -0.944   0.3486  
depth                -3.746e-04  6.284e-04  -0.596   0.5531  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01492 on 69 degrees of freedom
Multiple R-squared:  0.1063,    Adjusted R-squared:  0.04159 
F-statistic: 1.642 on 5 and 69 DF,  p-value: 0.1605

uomisto H (2010a). “A diversity of beta diver- sities: straightening up a concept gone awry. 1. Defining beta diversity as a function of alpha and gamma diversity.” Ecography, 33, 2–2

Prediction plots

Observed Species

predict(obsMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18       19 
219.0394 219.0394 191.6764 191.6764 200.2007 157.3011 157.3011 183.0517 202.9587 295.9379 295.9379 268.5750 268.5750 277.0992 234.1997 234.1997 259.9502 259.9502 327.4518 
      20       21       22       23       24       25       26       27       28       29       30       31       32       33       34       35       36       37       38 
327.4518 300.0888 300.0888 308.6131 265.7136 265.7136 291.4641 311.3712 311.3712 332.8183 332.8183 305.4553 305.4553 313.9796 313.9796 271.0800 271.0800 296.8306 296.8306 
      39       40       41       42       43       44       45       46       47       48       49       50       51       52       53       54       55       56       57 
321.6979 294.3349 294.3349 302.8592 302.8592 259.9597 259.9597 259.9597 285.7102 285.7102 305.6173 305.6173 290.3570 290.3570 262.9941 262.9941 271.5184 271.5184 228.6188 
      58       59       60       61       62       63       64       65       66       67       68       69       70       71       72       73       74       75 
228.6188 228.6188 254.3694 254.3694 274.2764 274.2764 249.3269 221.9639 221.9639 230.4882 230.4882 187.5886 187.5886 187.5886 213.3392 213.3392 233.2462 233.2462 

$se.fit
 [1] 26.41573 26.41573 25.69552 25.69552 25.55840 27.19115 27.19115 28.71794 33.31516 18.53028 18.53028 17.92699 17.92699 16.83554 19.61469 19.61469 21.04649 21.04649 18.84519
[20] 18.84519 18.39317 18.39317 16.87013 19.66040 19.66040 20.75488 27.57872 27.57872 18.84306 18.84306 18.36078 18.36078 16.60167 16.60167 19.19454 19.19454 20.13853 20.13853
[39] 18.16789 17.54242 17.54242 15.65854 15.65854 18.05154 18.05154 18.05154 19.02151 19.02151 25.99078 25.99078 18.90217 18.90217 18.03402 18.03402 16.32633 16.32633 18.03550
[58] 18.03550 18.03550 19.08494 19.08494 25.56051 25.56051 23.94049 23.00762 23.00762 21.87070 21.87070 22.64968 22.64968 22.64968 23.62789 23.62789 28.60636 28.60636

$df
[1] 69

$residual.scale
[1] 75.50836

diversityData$pred_obs = predict(obsMod, se.fit = TRUE)$fit
diversityData$se_obs = predict(obsMod, se.fit = TRUE)$se.fit

plotSpecs2 <- list(
  facet_wrap(~Depth, ncol = 1) ,
  theme_bw(base_size = 16) ,
  #geom_point(size = 4) ,
  geom_path(aes(color = as.factor(Station))) ,
  scale_x_log10(breaks = my_sizes, labels = as.character(my_sizes)) ,
  #scale_y_log10nice() ,
  scale_shape_manual(values = rep(21:25, 2)) ,
  scale_fill_viridis_d(option = "plasma") ,
  scale_color_viridis_d(option = "plasma") ,
  labs(x = expression(paste("Particle Size (", mu, "m)"))) ,
  theme(legend.position = "none",
        plot.margin = unit(c(0, 0, 0, 0), "cm"),
        axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = .5),
        axis.title.y = element_text(margin = unit(c(3, 3, 3, 3), "mm"), vjust = 0))
)

plotObs_pred <-  diversityData %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_obs, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_obs - 2 * se_obs, yend = pred_obs + 2 * se_obs, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted  ASVs") 
plotObs_pred

Richness

predict(chao1Mod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18       19 
414.2733 414.2733 319.9690 319.9690 396.7495 263.9882 263.9882 373.9864 352.6608 620.7464 620.7464 526.4421 526.4421 603.2225 470.4612 470.4612 580.4595 580.4595 703.5800 
      20       21       22       23       24       25       26       27       28       29       30       31       32       33       34       35       36       37       38 
703.5800 609.2758 609.2758 686.0562 553.2949 553.2949 663.2932 641.9675 641.9675 714.7463 714.7463 620.4421 620.4421 697.2225 697.2225 564.4612 564.4612 674.4595 674.4595 
      39       40       41       42       43       44       45       46       47       48       49       50       51       52       53       54       55       56       57 
681.7164 587.4121 587.4121 664.1925 664.1925 531.4312 531.4312 531.4312 641.4295 641.4295 620.1038 620.1038 592.5451 592.5451 498.2408 498.2408 575.0212 575.0212 442.2600 
      58       59       60       61       62       63       64       65       66       67       68       69       70       71       72       73       74       75 
442.2600 442.2600 552.2582 552.2582 530.9325 530.9325 477.2943 382.9901 382.9901 459.7705 459.7705 327.0092 327.0092 327.0092 437.0074 437.0074 415.6818 415.6818 

$se.fit
 [1]  86.97476  86.97476  84.60346  84.60346  84.15196  89.52786  89.52786  94.55488 109.69140  61.01161  61.01161  59.02528  59.02528  55.43162  64.58209  64.58209  69.29635
[18]  69.29635  62.04847  62.04847  60.56018  60.56018  55.54551  64.73259  64.73259  68.33621  90.80394  90.80394  62.04146  62.04146  60.45353  60.45353  54.66162  54.66162
[35]  63.19872  63.19872  66.30685  66.30685  59.81845  57.75905  57.75905  51.55631  51.55631  59.43537  59.43537  59.43537  62.62903  62.62903  85.57562  85.57562  62.23608
[52]  62.23608  59.37768  59.37768  53.75503  53.75503  59.38253  59.38253  59.38253  62.83786  62.83786  84.15892  84.15892  78.82495  75.75342  75.75342  72.01009  72.01009
[69]  74.57490  74.57490  74.57490  77.79571  77.79571  94.18749  94.18749

$df
[1] 69

$residual.scale
[1] 248.6141

diversityData$pred_chao1 = predict(chao1Mod, se.fit = TRUE)$fit
diversityData$se_chao1 = predict(chao1Mod, se.fit = TRUE)$se.fit

plotChao1_pred <-  diversityData %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_chao1, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_chao1 - 2 * se_chao1, yend = pred_chao1 + 2 * se_chao1, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"), alpha = 0.5)  +
  plotSpecs2 + ylab("Predictd Richness (Chao1)") + scale_y_log10()
plotChao1_pred

Shannon Diversity

predict(shanMod, se.fit = TRUE)

$fit
       1        2        3        4        5        6        7        8        9       10       11       12       13       14       15       16       17       18       19 
4.460000 4.460000 4.316870 4.316870 4.263512 4.006600 4.006600 4.096129 4.380305 4.922816 4.922816 4.779687 4.779687 4.726328 4.469417 4.469417 4.558945 4.558945 5.121512 
      20       21       22       23       24       25       26       27       28       29       30       31       32       33       34       35       36       37       38 
5.121512 4.978383 4.978383 4.925024 4.668113 4.668113 4.757641 5.041818 5.041818 5.170247 5.170247 5.027117 5.027117 4.973759 4.973759 4.716847 4.716847 4.806376 4.806376 
      39       40       41       42       43       44       45       46       47       48       49       50       51       52       53       54       55       56       57 
5.119395 4.976266 4.976266 4.922907 4.922907 4.665996 4.665996 4.665996 4.755524 4.755524 5.039701 5.039701 4.956219 4.956219 4.813089 4.813089 4.759731 4.759731 4.502820 
      58       59       60       61       62       63       64       65       66       67       68       69       70       71       72       73       74       75 
4.502820 4.502820 4.592348 4.592348 4.876525 4.876525 4.735050 4.591920 4.591920 4.538562 4.538562 4.281650 4.281650 4.281650 4.371179 4.371179 4.655356 4.655356 

$se.fit
 [1] 0.15808462 0.15808462 0.15377457 0.15377457 0.15295393 0.16272512 0.16272512 0.17186218 0.19937421 0.11089422 0.11089422 0.10728388 0.10728388 0.10075207 0.11738388
[16] 0.11738388 0.12595249 0.12595249 0.11277880 0.11277880 0.11007370 0.11007370 0.10095908 0.11765743 0.11765743 0.12420735 0.16504451 0.16504451 0.11276605 0.11276605
[31] 0.10987985 0.10987985 0.09935253 0.09935253 0.11486948 0.11486948 0.12051879 0.12051879 0.10872554 0.10498238 0.10498238 0.09370833 0.09370833 0.10802925 0.10802925
[46] 0.10802925 0.11383402 0.11383402 0.15554155 0.15554155 0.11311979 0.11311979 0.10792438 0.10792438 0.09770471 0.09770471 0.10793321 0.10793321 0.10793321 0.11421358
[61] 0.11421358 0.15296658 0.15296658 0.14327160 0.13768881 0.13768881 0.13088496 0.13088496 0.13554674 0.13554674 0.13554674 0.14140086 0.14140086 0.17119443 0.17119443

$df
[1] 69

$residual.scale
[1] 0.4518789

diversityData$pred_shanH = predict(shanMod, se.fit = TRUE)$fit
diversityData$se_shanH = predict(shanMod, se.fit = TRUE)$se.fit

plotShannonH_pred <- diversityData %>%

 #filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_shanH, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_shanH - 2 * se_shanH, yend = pred_shanH + 2 * se_shanH, xend = Size_Class, color = as.factor(Station)), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"),  alpha = 0.5)  +
  plotSpecs2 + ylab("Predicted Diversity (Shannon H)") #+ scale_y_log10()
plotShannonH_pred

Evenness

predict(pielouMod, se.fit = TRUE)

$fit
          1           2           3           4           5           6           7           8           9          10          11          12          13          14 
0.020303322 0.020303322 0.022872940 0.022872940 0.021919895 0.025509948 0.025509948 0.022656554 0.019388204 0.011062673 0.011062673 0.013632292 0.013632292 0.012679247 
         15          16          17          18          19          20          21          22          23          24          25          26          27          28 
0.016269300 0.016269300 0.013415905 0.013415905 0.006958323 0.006958323 0.009527941 0.009527941 0.008574897 0.012164950 0.012164950 0.009311555 0.006043205 0.006043205 
         29          30          31          32          33          34          35          36          37          38          39          40          41          42 
0.005735595 0.005735595 0.008305213 0.008305213 0.007352168 0.007352168 0.010942221 0.010942221 0.008088827 0.008088827 0.006506674 0.009076292 0.009076292 0.008123248 
         43          44          45          46          47          48          49          50          51          52          53          54          55          56 
0.008123248 0.011713300 0.011713300 0.011713300 0.008859906 0.008859906 0.005591556 0.005591556 0.009378051 0.009378051 0.011947669 0.011947669 0.010994624 0.010994624 
         57          58          59          60          61          62          63          64          65          66          67          68          69          70 
0.014584677 0.014584677 0.014584677 0.011731283 0.011731283 0.008462932 0.008462932 0.013402419 0.015972038 0.015972038 0.015018993 0.015018993 0.018609046 0.018609046 
         71          72          73          74          75 
0.018609046 0.015755652 0.015755652 0.012487301 0.012487301 

$se.fit
 [1] 0.005218045 0.005218045 0.005075779 0.005075779 0.005048692 0.005371218 0.005371218 0.005672814 0.006580929 0.003660388 0.003660388 0.003541218 0.003541218 0.003325617
[15] 0.003874598 0.003874598 0.004157430 0.004157430 0.003722594 0.003722594 0.003633304 0.003633304 0.003332450 0.003883628 0.003883628 0.004099827 0.005447777 0.005447777
[29] 0.003722173 0.003722173 0.003626906 0.003626906 0.003279421 0.003279421 0.003791603 0.003791603 0.003978075 0.003978075 0.003588804 0.003465251 0.003465251 0.003093117
[43] 0.003093117 0.003565821 0.003565821 0.003565821 0.003757425 0.003757425 0.005134104 0.005134104 0.003733850 0.003733850 0.003562360 0.003562360 0.003225030 0.003225030
[57] 0.003562651 0.003562651 0.003562651 0.003769953 0.003769953 0.005049109 0.005049109 0.004729098 0.004544822 0.004544822 0.004320241 0.004320241 0.004474116 0.004474116
[71] 0.004474116 0.004667349 0.004667349 0.005650773 0.005650773

$df
[1] 69

$residual.scale
[1] 0.01491559

diversityData$pred_pielouJ = predict(pielouMod, se.fit = TRUE)$fit
diversityData$se_pielouJ = predict(pielouMod, se.fit = TRUE)$se.fit

plot_pielouJ_pred <- diversityData %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ - 2 * se_pielouJ, yend = pred_pielouJ + 2 * se_pielouJ, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J)") + scale_y_log10()
plot_pielouJ_pred

Combined prediction plot

plotPredictions <- plot_grid(plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred, nrow = 1, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).

plotPredictions

ggsave(here::here("Figures", "ConventionalAlphaPredictions.png"), plotPredictions, width = 11, height = 4)

Even combindeder

plot_grid(plotObs, plotChao1, plotShan, plotPielou,
          plotObs_pred, plotChao1_pred, plotShannonH_pred, plot_pielouJ_pred,
          nrow = 2, labels = LETTERS)

Warning: NaNs producedWarning: Transformation introduced infinite values in continuous y-axisWarning: Removed 11 rows containing missing values (`geom_segment()`).

Combined summary table

alphaSummary <- tibble(
  metric = c("Observed ASVs", "Richness (Chao1)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(obsMod, chao1Mod, shanMod, pielouMod)
)

alphaSummary <- alphaSummary %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

alphaSummary <- alphaSummary %>%
  select(-model) %>%
  unnest(df)

alphaSummary <- alphaSummary %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()

alphaSummary %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>%
  bold(i = ~ p< 0.05, j = "p") %>%
  colformat_md() %>%
  set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")

Metric	Term	Estimate	Standard Error	T Value	p
Observed ASVs	Intercept	2.5 x 105	1.3 x 105	1.87	0.065
	log(Size Class)	3.4 x 101	8.0 x 100	4.19	< 0.001
	log(Size Class)2	-6.5 x 100	1.5 x 100	-4.33	< 0.001
	Latitude	-1.3 x 104	6.9 x 103	-1.87	0.066
	Latitude2	1.7 x 102	9.0 x 101	1.87	0.066
	Depth	3.9 x 100	3.2 x 100	1.23	0.224
Richness (Chao1)	Intercept	8.0 x 105	4.4 x 105	1.82	0.073
	log(Size Class)	9.0 x 101	2.6 x 101	3.40	0.001
	log(Size Class)2	-1.8 x 101	4.9 x 100	-3.61	< 0.001
	Latitude	-4.1 x 104	2.3 x 104	-1.82	0.073
	Latitude2	5.4 x 102	3.0 x 102	1.82	0.073
	Depth	1.7 x 101	1.0 x 101	1.62	0.110
Diversity (Shannon H)	Intercept	1.4 x 103	7.9 x 102	1.79	0.078
	log(Size Class)	2.1 x 10-1	4.8 x 10-2	4.28	< 0.001
	log(Size Class)2	-3.7 x 10-2	8.9 x 10-3	-4.14	< 0.001
	Latitude	-7.4 x 101	4.1 x 101	-1.78	0.079
	Latitude2	9.6 x 10-1	5.4 x 10-1	1.78	0.079
	Depth	1.3 x 10-2	1.9 x 10-2	0.69	0.490
Evenness (Pielou J)	Intercept	-2.5 x 101	2.6 x 101	-0.95	0.348
	log(Size Class)	-4.1 x 10-3	1.6 x 10-3	-2.61	0.011
	log(Size Class)2	7.1 x 10-4	3.0 x 10-4	2.40	0.019
	Latitude	1.3 x 100	1.4 x 100	0.94	0.348
	Latitude2	-1.7 x 10-2	1.8 x 10-2	-0.94	0.349
	Depth	-3.7 x 10-4	6.3 x 10-4	-0.60	0.553

Now considering breakaway values

richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.))))

Warning: `rename_()` was deprecated in dplyr 0.7.0.
Please use `rename()` instead.

diversityDataWB <- full_join(diversityData,
                             richSummary %>% rename_(.dots = setNames(names(.), paste0('break_', names(.)))),
                             by = c("ID" = "break_sample_names"), suffix = c("", "_break")) %>%
  mutate(pielouJ2 = shannonH/break_estimate) %>%
  identity()

diversityDataWB

pielouMod2 <- lm(pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + lat + I(lat^2) + depth, data = diversityDataWB)
summary(pielouMod2)

Call:
lm(formula = pielouJ2 ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.013935 -0.005053 -0.002494  0.000907  0.105945 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -1.980e+01  2.751e+01  -0.720   0.4743  
log(Size_Class)      -3.292e-03  1.667e-03  -1.975   0.0523 .
I(log(Size_Class)^2)  5.754e-04  3.102e-04   1.855   0.0679 .
lat                   1.030e+00  1.434e+00   0.718   0.4751  
I(lat^2)             -1.338e-02  1.866e-02  -0.717   0.4760  
depth                -2.432e-04  6.603e-04  -0.368   0.7138  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01567 on 69 degrees of freedom
Multiple R-squared:  0.06742,   Adjusted R-squared:  -0.0001586 
F-statistic: 0.9977 on 5 and 69 DF,  p-value: 0.4258

Ok. So the narrative makes sense. Alpha diveristy is driven by variability in richness rather than evenness. Why would we see an effect in chao1 but not breakaway? Because chao1 is more driven by abundant stuff that makes the rarification threshold. My first hunch is that chao1 responds to evenness, but actually that shouldn’t have any effect since there is no evenness variability? Or maybe just that breakaway is more variable (because it detects fine level differences in rare species) and that doesn’t map as nicely with overall patterns.

plotBreak <- diversityDataWB %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Richness (Breakaway)")
plotBreak

plotPielou2 <- diversityDataWB %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  #scale_y_log10()+
  ylab("Evenness (PielouJ)")
plotPielou2

Redo predictions for good measure

predict(pielouMod2, se.fit = TRUE)

$fit
            1             2             3             4             5             6             7             8             9            10            11            12 
 0.0116609352  0.0116609352  0.0135814780  0.0135814780  0.0133501419  0.0160295241  0.0160295241  0.0139649206  0.0099302928  0.0042917835  0.0042917835  0.0062123262 
           13            14            15            16            17            18            19            20            21            22            23            24 
 0.0062123262  0.0059809902  0.0086603723  0.0086603723  0.0065957688  0.0065957688  0.0010653312  0.0010653312  0.0029858739  0.0029858739  0.0027545378  0.0054339200 
           25            26            27            28            29            30            31            32            33            34            35            36 
 0.0054339200  0.0033693165 -0.0006653113 -0.0006653113  0.0001751617  0.0001751617  0.0020957044  0.0020957044  0.0018643683  0.0018643683  0.0045437505  0.0045437505 
           37            38            39            40            41            42            43            44            45            46            47            48 
 0.0024791470  0.0024791470  0.0008731389  0.0027936816  0.0027936816  0.0025623455  0.0025623455  0.0052417277  0.0052417277  0.0052417277  0.0031771242  0.0031771242 
           49            50            51            52            53            54            55            56            57            58            59            60 
-0.0008575036 -0.0008575036  0.0032944705  0.0032944705  0.0052150133  0.0052150133  0.0049836772  0.0049836772  0.0076630594  0.0076630594  0.0076630594  0.0055984559 
           61            62            63            64            65            66            67            68            69            70            71            72 
 0.0055984559  0.0015638281  0.0015638281  0.0066369468  0.0085574895  0.0085574895  0.0083261534  0.0083261534  0.0110055356  0.0110055356  0.0110055356  0.0089409321 
           73            74            75 
 0.0089409321  0.0049063043  0.0049063043 

$se.fit
 [1] 0.005482717 0.005482717 0.005333235 0.005333235 0.005304774 0.005643660 0.005643660 0.005960553 0.006914730 0.003846052 0.003846052 0.003720837 0.003720837 0.003494300
[15] 0.004071127 0.004071127 0.004368305 0.004368305 0.003911413 0.003911413 0.003817594 0.003817594 0.003501480 0.004080615 0.004080615 0.004307780 0.005724101 0.005724101
[29] 0.003910971 0.003910971 0.003810871 0.003810871 0.003445761 0.003445761 0.003983922 0.003983922 0.004179853 0.004179853 0.003770837 0.003641017 0.003641017 0.003250008
[43] 0.003250008 0.003746688 0.003746688 0.003746688 0.003948011 0.003948011 0.005394518 0.005394518 0.003923240 0.003923240 0.003743051 0.003743051 0.003388611 0.003388611
[57] 0.003743358 0.003743358 0.003743358 0.003961175 0.003961175 0.005305212 0.005305212 0.004968969 0.004775346 0.004775346 0.004539374 0.004539374 0.004701055 0.004701055
[71] 0.004701055 0.004904088 0.004904088 0.005937394 0.005937394

$df
[1] 69

$residual.scale
[1] 0.01567214

diversityDataWB$pred_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$fit
diversityDataWB$se_pielouJ2 = predict(pielouMod2, se.fit = TRUE)$se.fit

plot_pielouJ2_pred <- diversityDataWB %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = pred_pielouJ2, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_pielouJ2 - 2 * se_pielouJ2, yend = pred_pielouJ2 + 2 * se_pielouJ2, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Evenness (Pielou J2)") #+ scale_y_log10()
plot_pielouJ2_pred

Breakaway richness subplots

plotBreakaway <- diversityDataWB %>%
#filter(Depth %in% c("Surface", "Bottom")) %>%
  ggplot(aes(x = Size_Class, y = break_estimate, shape = as.factor(Station), fill = as.factor(Station))) +
  plotSpecs +
  geom_errorbar(aes(ymin = break_lower, ymax = break_upper), width = -.1) + 
  scale_y_log10() +
  ylab("Richness (Breakaway)")
plotBreakaway

#predict(breakLm, se.fit = TRUE)
# doesn't work because built with a different data frame

Why are these not smooth curves?!! What if I redo the model, this time with the same data frame

breakLm2 <- lm(break_estimate ~ log(Size_Class) + I(log(Size_Class) ^2) + lat +  I(lat^2) + depth ,data = diversityDataWB)
breakLm2 %>% summary()

Call:
lm(formula = break_estimate ~ log(Size_Class) + I(log(Size_Class)^2) + 
    lat + I(lat^2) + depth, data = diversityDataWB)

Residuals:
    Min      1Q  Median      3Q     Max 
-2974.5 -1191.2  -151.6   599.9  6768.1 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)  
(Intercept)          7124615.61 3339862.88   2.133   0.0365 *
log(Size_Class)          244.45     202.35   1.208   0.2312  
I(log(Size_Class)^2)     -75.16      37.65  -1.996   0.0498 *
lat                  -370568.38  174035.93  -2.129   0.0368 *
I(lat^2)                4817.28    2265.61   2.126   0.0371 *
depth                    151.10      80.15   1.885   0.0636 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1902 on 69 degrees of freedom
Multiple R-squared:  0.1414,    Adjusted R-squared:  0.0792 
F-statistic: 2.273 on 5 and 69 DF,  p-value: 0.0567

Note the non statistical significance overall

#predict(breakLm2, se.fit = TRUE)
diversityDataWB$pred_break = predict(breakLm2, se.fit = TRUE)$fit
diversityDataWB$se_break = predict(breakLm2, se.fit = TRUE)$se.fit

plot_break_pred <- diversityDataWB %>%

#filter(Depth %in% c("Surface", "Bottom")) %>%
#  filter(Station == 4.3) %>%
  ggplot(aes(x = Size_Class, y = pred_break, shape = as.factor(Station), fill = as.factor(Station))) +
  geom_segment(aes(y = pred_break - 2 * se_break, yend = pred_break + 2 * se_break, xend = Size_Class, color = as.factor(Station), alpha = 0.5), arrow = arrow(angle = 70, length = unit(0.05, "in"), ends = "both"))  +
  plotSpecs2 + ylab("Predicted Richness (Breakaway -- LM)") #+ scale_y_log10()
plot_break_pred

Rebuilding combined products

plotAlphaWB <- plot_grid(plotBreakaway, plotShan, plotPielou2, nrow = 1, labels = LETTERS)
plotAlphaWB

ggsave(here::here("Figures", "BreakawayAlpha.png"), plotAlpha, width = 11, height = 4)

Summary table I want both breakaway metrics here

bettaTable <- myBet$table %>% 
  as.data.frame() %>%
  rename(estimate = Estimates,
         `std.error` = `Standard Errors`,
         `p.value`=`p-values`
         ) %>%
  mutate(`statistic` = NA) %>%
  rownames_to_column(var = "term") %>%
  select(term, estimate, std.error, statistic, p.value) %>%
  as_tibble()
bettaTable

alphaSummary2 <- tibble(
  metric = c("Richness (Breakaway -- LM)", "Diversity (Shannon H)", "Evenness (Pielou J)"),
  model = list(breakLm, shanMod, pielouMod2)
)
  
alphaSummary2 <- alphaSummary2 %>%
  mutate(df = map(model, ~broom::tidy(summary(.))))

## Add in willis variables

breakawaySummary <- tibble(
  metric = "Richness (Breakaway -- Betta)",
  model = NULL,
  df = list(bettaTable)
)

alphaSummary2 = bind_rows(breakawaySummary, alphaSummary2)

alphaSummary2 <- alphaSummary2 %>%
  select(-model) %>%
  unnest(df)

alphaSummary2 <- alphaSummary2 %>%
  rename(Metric = metric, Term = term, Estimate = estimate, `Standard Error` = std.error, `T Value` = statistic, p = p.value) %>%
  mutate(Term = str_replace(Term, "^I?\\((.*)\\)", "\\1"),
         Term = str_replace(Term, "\\^2", "\\^2\\^"),
         Term = str_replace(Term, "depth", "Depth"),
         Term = str_replace(Term, "lat", "Latitude"),
         Term = str_replace(Term, "_", " ")# BOOKMARK!!
         ) %>%
  mutate(Estimate = format(Estimate, digits = 2, scientific = TRUE) %>%
           reformat_sci()
         ) %>%
  mutate(`Standard Error` = format(`Standard Error`, digits = 2, scientific = TRUE) %>%
           reformat_sci()
  ) %>%
  mutate(`T Value` = format(`T Value`, digits = 2, scientific = FALSE)) %>%
  mutate(p = if_else(p < 0.001, "< 0.001", format(round(p, digits = 3)))) %>%
  rename(`Standard\nError` = `Standard Error`) %>%
  identity()



alphaSummary2

alphaTable2 <- alphaSummary2 %>% flextable() %>% merge_v(j = 1) %>% theme_vanilla() %>% bold(i = ~ p< 0.05, j = "p") %>% colformat_md() %>% set_table_properties(layout = "autofit") %>%
  align(j = -c(1:2), align = "right")
alphaTable2

Metric	Term	Estimate	Standard Error	T Value	p
Richness (Breakaway – Betta)	Intercept	7.1 x 106	2.4 x 102	NA	< 0.001
	log(Size Class)	1.2 x 102	6.1 x 101	NA	0.058
	log(Size Class)2	-5.0 x 101	1.2 x 101	NA	< 0.001
	Latitude	-3.7 x 105	6.1 x 100	NA	< 0.001
	Latitude2	4.8 x 103	1.6 x 10-1	NA	< 0.001
	Depth	1.5 x 102	1.0 x 101	NA	< 0.001
Richness (Breakaway – LM)	Intercept	7.1 x 106	3.3 x 106	2.13	0.036
	log(Size Class)	2.4 x 102	2.0 x 102	1.21	0.231
	log(Size Class)2	-7.5 x 101	3.8 x 101	-2.00	0.050
	Latitude	-3.7 x 105	1.7 x 105	-2.13	0.037
	Latitude2	4.8 x 103	2.3 x 103	2.13	0.037
	Depth	1.5 x 102	8.0 x 101	1.89	0.064
Diversity (Shannon H)	Intercept	1.4 x 103	7.9 x 102	1.79	0.078
	log(Size Class)	2.1 x 10-1	4.8 x 10-2	4.28	< 0.001
	log(Size Class)2	-3.7 x 10-2	8.9 x 10-3	-4.14	< 0.001
	Latitude	-7.4 x 101	4.1 x 101	-1.78	0.079
	Latitude2	9.6 x 10-1	5.4 x 10-1	1.78	0.079
	Depth	1.3 x 10-2	1.9 x 10-2	0.69	0.490
Evenness (Pielou J)	Intercept	-2.0 x 101	2.8 x 101	-0.72	0.474
	log(Size Class)	-3.3 x 10-3	1.7 x 10-3	-1.97	0.052
	log(Size Class)2	5.8 x 10-4	3.1 x 10-4	1.86	0.068
	Latitude	1.0 x 100	1.4 x 100	0.72	0.475
	Latitude2	-1.3 x 10-2	1.9 x 10-2	-0.72	0.476
	Depth	-2.4 x 10-4	6.6 x 10-4	-0.37	0.714

alphaTable2 %>% save_as_docx(path = here::here("Tables", "alphaTable2.docx"))

myBet$table

And finally predictions from richness, diversity evenness again.

plot_grid(plot_break_pred,plotShannonH_pred,plot_pielouJ2_pred, nrow = 1, labels = LETTERS)